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ABSTRACT 

HEALPix - the Hierarchical Equal Area iso-Latitude Pixelization - is a ver- 
satile data structure with an associated library of computational algorithms and 
visualization software that supports fast scientific applications executable directly 
on very large volumes of astronomical data and large area surveys in the form of 
discretized spherical maps. Originally developed to address the data processing 
and analysis needs of the present generation of cosmic microwave background 
(CMB) experiments (e.g. BOOMERanG, WMAP), HEALPix can be expanded 
to meet many of the profound challenges that will arise in confrontation with the 
observational output of future missions and experiments, including e.g. Planck, 
Herschel, SAFIR, and the Beyond Einstein CMB polarization probe. In this pa- 
per we consider the requirements and constraints to be met in order to implement 
a sufficient framework for the efficient discretization and fast analysis/synthesis 
of functions defined on the sphere, and summarise how they are satisfied by 
HEALPix. 

Subject headings: cosmic microwave background — cosmology: observations - 
methods: statistical 
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1. Introduction 



Advanced detectors in modern astronomy produce data at huge rates at many wave- 
lengths. Some data sets are becoming very impressively large indeed. Of particular interest 
to us are those that accumulate astronomical data distributed on the entire sky or a con- 
siderable fraction thereof. Typical examples include radio, cosmic microwave background, 
submillimeter, infrared, X-ray, and gamma-ray sky maps of diffuse emission, and full sky or 
wide area surveys of extragalactic objects. Together with the wealth of gathered informa- 
tion comes the inevitable burden of increased complexity of the tasks of data reduction and 
science extraction. In this paper we are focused on those issues, which are related to the 
distinctive nature of the spherical spatial domain over which the data reside, and need to be 
analysed for scientific return. Our original motivations come from work in the field of mea- 
surement and interpretation of the cosmic microwave background (CMB) anisotropy. The 
growing complexity of the CMB anisotropy science extraction problem can be illustrated by 
the transition between the data sets of COBE-DMR (early 1990s, 7 deg resolution, 6000 pixel 
sky maps at 3 wavelengths), Boomerang (late 1990s, FWHM 12 arcmin, partial sky maps of 
200000 pixels at 4 wavelengths), WMAP (early 2000s, resolution up to FWHM 14 arcmin, 
3 million pixel sky maps at 5 wavelengths), and Planck (data expected 2008, resolution up 
to FWHM 5 arcmin, 50 million pixel sky maps at 9 wavelengths). Science extraction from 
these data sets involves (1) global analysis problems - harmonic decomposition, estimation 
of the power spectrum and higher order measures of spatial correlations, (2) simulations of 
models of the primary and foreground sky signals to study instrument performance, and 
calibrate foreground separation and statistical inference methods, (3) real space morpholog- 
ical analyses, object detection, identification, and characterization, and (4) spatial and/or 
spectral cross-correlation with external data sets. These tasks, and many others, necessitate 
a careful definition of the data models, and proper set-up of the mathematical framework of 
data analysis such that the algorithmic and computing time requirements can be satisfied in 
order to achieve successful scientific interpretation of expensive and precious observations. 
A particular method of addressing some of these problems is described next. 



2. Discretized Mapping and Analysis of Functions on the Sphere 

The analysis of functions on domains with spherical topology occupies a central place 
in physical science and engineering disciplines. This is particularly apparent in the fields of 
astronomy, cosmology, geophysics, atomic and nuclear physics. In many cases the geometry 
is either dictated by the object under study or approximate spherical symmetry can be 
exploited to yield powerful perturbation methods. Practical limits for the purely analytical 
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study of these problems create an urgent necessity for efficient and accurate numerical tools. 

The simplicity of the spherical form belies the intricacy of global analysis on the sphere. 
There is no known point set which achieves the analogue of uniform sampling in Euclidean 
space and allows exact and invertible discrete spherical harmonic decompositions of arbitrary 
but band-limited functions. Any existing proposition of practical schemes for the discrete 
treatment of such functions on the sphere introduces some (hopefully small) systematic error 
dependent on the global properties of the point set. The goal is to minimise these errors and 
faithfully represent deterministic functions as well as realizations of random variates both in 
configuration and Fourier space while maintaining computational efficiency 

We illustrate these points using as an example the field of Cosmic Microwave Background 
(CMB) anisotropies. Here we are already in a situation of an ongoing rapid growth of the 
volume of the available data. The NASA's Wilkinson Microwave Anisotropy Probe (WMAP) 
already is, and ESA's Planck will (in the near future) be aiming to provide multi-frequency, 
high resolution, full sky measurements of the anisotropy in both temperature and polarization 
of the cosmic microwave background radiation. The ultimate data products of these missions 
- multiple microwave sky maps, each of which will have to comprise more than ~ 10 6 pixels 
in order to render the angular resolution of the instruments — will present serious challenges 
to those involved in the analysis and scientific exploitation of the results of both surveys. 

A digitized sky map is an essential intermediate stage in information processing between 
the entry point of data acquisition by the instruments - - very large time ordered data 
streams, and the final stage of astrophysical analysis — typically producing a 'few' numerical 
values of physical parameters of interest. COBE-DMR sky maps (angular resolution of 7° 
(FWHM) in three frequency bands, two channels each, 6144 pixels per map) were considered 
large at the time of their release. As for the presently available (WMAP) and future CMB 
maps, a whole sky CMB survey at the angular resolution of ~ 10' (FWHM), discretized 
with a few pixels per resolution element (so that the discretisation effects on the signal are 
sub-dominant with respect to the effects of instrument's angular response), will require map 
sizes of at least N pix ~ a few xl.5 10 6 pixels. Even more pixels than that will be needed to 
represent the Planck-HFI higher resolution high-frequency channels. 

This estimate, N pix , should be multiplied by the number of frequency bands (or, indeed, 
by the number of individual observing channels — 74 in the case of Planck — for the analysis 
work to be done before the final coadded maps are made for each frequency band) to render 
an approximate expected size of the already very compressed form of survey data which 
would be the input to the astrophysical analysis pipeline. 

Hence, a careful attention ought to be given to devising high resolution CMB map 
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structures which can maximally facilitate the forthcoming analyses of large size data sets, 
especially so because many essential scientific questions can only be answered by global 
studies of such data sets. 

This paper describes the essential geometric and algebraic properties of our method 
of digital representation of functions on the sphere - - the Hierarchical Equal Area and 
isoLatitude Pixelization (HEALPix) — and the associated multi-purpose computer software 
package. We have originally devised HEALPix , and we started distributing HEALPix 
software to the community in 1997. Presently HEALPix software is distributed from the 
web-site www.eso.org/ science/healpix. 

3. Requirements for a Spherical Pixelization Scheme 

Numerical analysis of functions on the sphere involves (1) a class of mathematical oper- 
ations, whose objects are (2) discretised maps, i.e. quantizations of arbitrary functions ac- 
cording to a chosen tessellation (exhaustive partition of the sphere into finite area elements). 
Hereafter we mostly specialise our discussion to CMB related applications of HEALPix , but 
all our statements hold true generally for any relevant deterministic and random functions 
on the sphere. 

Considering point (1): Standard operations of numerical analysis which one might wish 
to execute on the sphere include convolutions with local and global kernels, Fourier analysis 
with spherical harmonics and power spectrum estimation, wavelet decomposition, nearest- 
neighbour searches, topological analysis, including searches for extrema or zero-crossings, 
computing Minkowski functionals, extraction of patches and finite differencing for solving 
partial differential equations. Some of these operations become prohibitively slow if the 
sampling of functions on the sphere, and the related structure of the discrete data set, are 
not designed carefully. 

Regarding point (2): Typically, a whole sky map rendered by a CMB experiment con- 
tains (i) signals coming from the sky, which are by design strongly band-width limited (in the 
sense of spatial Fourier decomposition) by the instrument's angular response function, and 
(ii) a projection into the elements of a discrete map, or pixels, of the observing instrument's 
noise; this pixel noise should be random, and white, at least near the discretisation scale, 
with a band-width significantly exceeding that of all the signals. 

With these considerations in mind we proposed the following list of desiderata for the 
mathematical structure of discrete whole sky maps: 
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1. Hierarchical structure of the database. This is recognized as essential for very 
large data bases, and was indeed postulated already in construction of the Quadrilat- 
eralized Spherical Cube (or QuadCube, see White & Stemwedel (1992), and 
http://lambda.gsfc.nasa.gov/product/cobe/skymap_info_new.cfm), which was used for 
all the COBE data. A simple argument in favour of this states that the data elements 
which are nearby in a multi-dimensional configuration space (here, on the surface of 
a sphere), are also nearby in the tree structure of the data base. This property fa- 
cilitates various topological methods of analysis, and allows for easy construction of 
wavelet transforms on triangular and quadrilateral grids through fast look-up of nearest 
neighbors. 

2. Equal areas of discrete elements of partition. This is advantageous because 
white noise at the sampling frequency of the instrument gets integrated exactly into 
white noise in the pixel space, and sky signals are sampled without regional dependence 
(although still care must be taken to choose a pixel size sufficiently small compared 
to the instrumental resolution to avoid excessive, and pixel shape dependent, signal 
smoothing). 

3. Iso-Latitude distribution of discrete area elements on a sphere. This property 
is essential for computational speed in all operations involving evaluations of spher- 
ical harmonics. Since the associated Legendre polynomials are evaluated via slow 
recursions, any sampling grid deviations from an iso-latitude distribution result in a 
prohibitive loss of computational performance with the growing number of sampling 
points. 

Various known sampling distributions on a sphere have been used for the discretisation 
and analysis of functions (for example, see Driscoll & Healy (1994), Muciaccia, Natoli & 
Vittorio (1998) — rectangular grids, Baumgardner & Frederickson (1985), Tegmark (1996) 
- icosahedral grids, Saff & Kuijlaars (1997), Crittenden & Turok (1998) — 'igloo' grids, 
and Szalay & Brunner (1998) — a triangular grid), but each fails to meet simultaneously all 
of these requirements. In particular, 

i) Quadrilateralized Spherical Cube obeys points 1 and (approximately) 2, but fails on 
point 3, and cannot be used for efficient Fourier analysis at high resolution. 

ii) Equidistant Cylindrical Projection, a very common computational tool in geophysics, 
and climate modeling, and recently implemented for CMB work (Muciaccia, Natoli, Vittorio, 
1998), satisfies points 1 and 3, but by construction fails with point 2. This is a nuisance 
from the point of view of application to full sky survey data due to wasteful oversampling 
near the poles of the map (or - while angular resolution of the measurements is fixed by the 
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instrument and does not vary over the sky, the map resolution, or pixel size, depends on the 
distance form the poles). 

iii) Hexagonal sampling grids with icosahedral symmetry perform superbly in those 
applications where near uniformity of sampling on a sphere is essential (Saff, Kuijlaars, 
1997), and can be devised to meet requirement 2 (see e.g. Tegmark 1996). However, by 
construction they fail both points 1 and 3. 

iv) Igloo- type constructions are devised to satisfy point 3 (E. Wright, 1997, private 
communication; Crittenden & Turok, 1998). Point 2 can be satisfied to reasonable accuracy 
if quite a large number of base-resolution pixels is used( which, however, precludes efficient 
construction of simple wavelet transforms). Conversely, a tree-structure seeded with a small 
number of base-resolution pixels forces significant variations in both area and shape of the 
pixels. 

v) GLESP construction (Doroshkevich etal 2003) takes advantage of the Gauss- Legendre 
quadrature to render high accuracy of numerical integration, but allows irregular variations 
the in pixel area and is not hierarchical — in fact it offers no relation between the tessellations 
derived at different resolutions. 



4. Meeting the Requirements 

All the requirements introduced in section 3. are satisfied by the class of spherical 
tessellations structured as follows. 

First let us assume that the sphere is partitioned into a number of curvilinear quadrilat- 
erals, which constitute the base level tessellation. If there exists a mapping of each element of 
partition onto a square [0, 1] x [0, 1], a nested nxn subdivision of the square into ever dimin- 
ishing sub-elements obtains trivially, and a hierarchical tree structure of resulting database 
follows. For example, a 2 x 2 partition renders a quadrilateral tree, which admits an el- 
egant binary indexation (illustrated in Fig. 1) previously employed in construction of the 
QuadCube spherical pixelization. 

Next, let us consider the base level spherical tessellation. An entire class of such tes- 
sellations can be constructed as illustrated in Figs. 2, and 3. These constructions are 
characterized by two parameters: Ng - a number of the base-resolution pixel layers between 
the north and south poles, and - a multiplicity of the meridional cuts, or the number 
of equatorial, or circum-polar base-resolutions pixels. Obviously, the total number of base- 
resolution pixels is equal to N pix = N 9 x N^, and the area of each one of them is equal to 
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Qpix = ^n/Ng/Ntj,. One may also notice (see Fig. 2) that each tessellation includes two single 
layers of polar cup pixels (with or without an azimuthal twist in their respective positions on 
the sphere for odd, and even values of Ng, respectively), and an Ng — 2 layers of equatorial 
zone pixels, which form a regular rhomboidal grid in the cylindrical projection of the sphere. 
Since the cylindrical projection is an area preserving mapping, this property immediately 
illustrates that the areas of equatorial zone pixels are all equal, and to meet our requirement 
of fully equal area partition of the sphere, we need to demonstrate that our constructions 
render identical areas of the polar pixels as well. Indeed, this allows to formulate a constraint 
on the latitude 9* at which the lateral vertices of both polar and equatorial pixels meet: 

2vr x (1 - cos0*) = N^x Q pix /2, hence cos 6^ = (N g - 1)/N g . (1) 

The curvilinear quadrilateral pixels of this tessellation class retain equal areas, but vary 
in shape depending on their positions on the sphere. We have chosen the Ng = 3, and = 4 
grid (middle row, right column in both Figs. 2, and 3) for definition of our digital full sky 
map data standard, and named it the HEALPix grid (Gorski et al. 1998). This grid is the 
basis for development of the HEALPix software, and it is described in more detail in the 
next section. 



5. The HEALPix Grid 

The requirements formulated in Sec. 3 are satisfied by construction with the Hierarchical 
Equal Area, iso-Latitude Pixelisation (HEALPix) of the sphere, which is shown in Figure 
(4). This Figure demonstrates that HEALPix is built geometrically as a self-similar, refinable 
quadrilateral mesh on the sphere. The base-resolution comprises twelve pixels in three 
rings around the poles and equator. The resolution of the grid is expressed by parameter 
N si( ie which defines the number of divisions along the side of a base-resolution pixel that 
is needed to reach a desired high-resolution partition. All pixel centers are placed on rings 
of constant latitude, and are equidistant in azimuth (on each ring). All iso-latitude rings 
located between the upper and lower corners of the equatorial base-resolution pixels (i.e. 
— 2/3 < cos# + < 2/3), or in the equatorial zone, are divided into the same number of pixels: 
N eq = AxN side . The remaining rings are located within the polar cap regions (| cos0*| > 2/3) 
and contain a varying number of pixels, increasing from ring to ring, with increasing distance 
from the poles, by one pixel within each quadrant. 

A HEALPix map has A pix = 12A s 2 ide pixels of the same area f2 pix = ^ . 
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5.1. Pixel Positions 

Locations on the sphere are defined by (z = cos#, 0) where 9 G [0, ir] is the colatitude in 
radians measured from the North Pole and G [0, 2ir] is the longitude in radians measured 
Eastward. 

For a resolution parameter A^e, the pixels are layed out on 4A / S id e — 1 iso-latitude rings, 
and can be ordered using the pixel index p G [0, N P i x ] running around those rings from the 
North to South pole. 

Pixel centers on the northern hemisphere are given by the following equations: 

North polar cap — for ph = (p+l)/2, the ring index 1 < i < N s i de , and the pixel-in-ring 
index 1 < j < 4i 



* = I{\Jp h -^Ph)) + l (2) 
3 = p+1-2^-1) (3) 



i 2 



Z = (4) 
= J0'-|)> and s = l (5) 

North equatorial belt — for p' = p — 2N sidc (N S i de — 1), iV S i d e < i < 2N side , and 1 < j < 
4Ar s idc 

i = I(p'/AN sidc ) + N sidc (6) 

3 = (p' mod 4iV sidc ) + 1 (7) 
4 2i 
3 SA'sidc 

= 77^- -Pj, and s = (i - N sidc + 1) mod 2, (9) 

^ v sidc v ^ 7 

where the auxilliary index s describes the phase shifts along the rings, and l(x) is the largest 
integer number smaller than x. 

Pixel center positions on the Southern hemisphere are obtained by the mirror symmetry 
of the grid with respect to the Equator (z — 0). 

Defining Az and A0 as the variation of z and when i and j are respectively increased 
by unity, one can check that discretized area element |AzA0| = f2 P i x , i.e. it is a constant. 
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5.2. Pixel Indexing 

Specific geometrical properties allow HEALPix to support two different numbering 
schemes for the pixels, as illustrated in the Figure (5). 

First, one can simply count the pixels moving down from the north to the south pole 
along each iso-latitude ring. It is in this scheme that Fourier transforms with spherical 
harmonics are easy to implement. 

Second, one can replicate the tree structure of pixel numbering used e.g. with the 
Quadrilateralized Spherical Cube. This can easily be implemented since, due to the simple 
description of pixel boundaries, the analytical mapping of the HEALPix base-resolution 
elements (curvilinear quadrilaterals) into a [0,1] x [0,1] square exists. This tree structure, 
a.k.a. nested scheme, allows one to implement efficiently all applications involving nearest- 
neighbour searches (see Wandelt, Hivon, and Gorski 1998), and also allows for an immediate 
construction of the fast Haar wavelet transform on HEALPix . 

The base-resolution pixel index number / runs in {0, NgN^ — 1} = {0, 11}. Introducing 
the row index 

/row = i(//i\g (10) 

we define two functions which index the location of the southernmost corner (or vertex) of 
each base resolution pixel on the sphere, in latitude and longitude respectively: 

Fl(f) = /row + 2, (11) 

F 2 (f) = 2(/ mod N+) - (/ row mod 2) + 1. (12) 

Consider the nested index p n e [0, 12-/V s 2 ide — 1], and define^ = (p n mod A" s 2 ide ), where p' n 
denotes the nested pixel index within each base-resolution element., which has the following 
binary representation p' n — [... (and bi = 0, or 1, and has the weight 2 l ). 

Given a grid resolution parameter A^de the location of a pixel center on each base- 
resolution pixel is represented by the 2 indices x and y e {0, A^dc — 1}- They both have 
their origin in the southernmost corner of each base-resolution pixel, with the x index running 
along the North-East direction, while the y index runs along the North- West direction. The 
binary representation of p' n determines the values of x, and y as the following combinations 
of even and odd bits, respectively 



x = [...b 2 b ] 2 
y = [...& 3 &i] 2 



(13) 
(14) 



-10- 



Next, we introduce the vertical and horizontal coordinates (within the base- resolution pixel) 

v = x + y (15) 
h = x — y (16) 

and obtain the following relations for the ring index % (e {1, (N + l)iV sidc — 1}) 

i = F 1 (f)N side -v-l, (17) 

and the longitude index 

. = F 2 (f)N sidc + h + s (ig) 
which can be translated into (z, <fi) coordinates using Eqs. 4, 5, 8, and 9. 

5.3. Pixel Boundaries 

Pixel boundaries are non-geodesic and take a very simple form: cos 6 = a + b x (f) in 
the equatorial zone, and cos# = a + b/<f) 2 in the polar caps. This allows one to explicitly 
check by simple analytical integration the exact area equality among pixels, and to compute 
efficiently more complex objects, e.g. the Fourier transforms of individual pixels. 

Since the pixel center location is parametrized by integer value of j, setting j — k + 1/2 
orj = k + 1/2 + i with k a positive integer in Eq. (5) and substituting into Eq. (4) will give 
for the pixels boundaries of the North polar cap (z > 2/3) 

^ = l-^fTTT^V. (20) 



where 4>t = <f) mod tt/2. The base pixels have boundaries defined as 



z>\, ( t ) = k % (21) 



with k' = 0,1,2,3 



Similarly, for 2/3 < z < —2/3 the pixel boundaries can be found by setting j = 
k + s/2 ±(i- N sidc )/2 in Eq. (9) and substituting into Eq. (8) 

~ 3 3iV sidc ± Stt 
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Using these pixel boundaries, one can easily check by integration that each individual 
pixel has the same surface area fi P i x - 

Table 1 summaries the number of pixels and resolution available in HEALPix . Since 
all pixels have the same surface area but slightly different shape, the angular resolution is 
defined as 

(24 ) 



6. Spherical Harmonic Transforms 

The requirement of iso-latitudinal location of pixel centers was built into HEALPix in 
order for the grid to support the fast discrete spherical harmonic transforms. The reason for 
the fast (Np{x) scaling of the harmonic transform computational time with the size of the grid 
(or the grid resolution) is entirely geometrical - the associated Legendre function components 
of spherical harmonics, which can only be generated via slow recursions, have to be evaluated 
only once for each pixel ring. For other grids, which are not iso-latitude constrainted, the 
extra computing time wasted on non-optimal generation of the associated Legendre functions, 
typically results with computational performance of order ~ N^ ix . This geometrical aspect 
of discrete spherical transform computation is illustrated in Fig. 6, which compares HEALPix 
with other tessellations including Quadrilateralized Spherical Cube (or QuadCube, used by 
NASA as data structure for COBE mission products), icosahedral tessellation of the sphere, 
and Equidistant Coordinate Partition, or the "geographic grid." This plot makes it visually 
clear why the iso-latitude ECP and HEALPix points-sets support faster computation of 
spherical harmonic transforms than the QuadCube, the icoshedral grid, and any non iso- 
latitude construction. 

Figure (7) demonstrates the fundamental difference between computing speeds, which 
can be achieved on iso-latitude and non iso-latitude point-sets. In order to be able to perform 
the necessary computational work in support of the multi-million pixel spherical data sets 
one has to resort to iso-latitude structures of point-sets/sky maps, e.g. HEALPix. Moreover, 
the future needs are already fairly clear - measurements of the CMB polarization will require 
massively multi-element arrays of detectors, and will produce data sets characterized by a 
great multiplicity (of order of a few thousand) of sky maps. Since there are no computation- 
ally faster methods than those already employed in HEALPix, and global synthesis/analysis 
of a multimillion pixel map consumes about 103s of a standard serial machine CPU time, 
the necessary speed-up will have to be achieved via optimized parallelization of the required 
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computing. 

A detailed description of implementation of the spherical harmonic transforms in the 
HEALPix software package, and the analysis of performance and accuracy thereof will given 
in a dedicated publication. 

7. Summary 

The Hierarchical Equal Area iso-Latitude Pixelization, HEALPix, is a methodology of 
discretization and fast numerical analysis and synthesis of functions or data distributed on 
the sphere. HEALPix is an intermediate data-structural, algorithmic, and functional layer 
between astronomical data, and the domain of application of variety of science extraction 
tools. HEALPix as a sky map format and associated set of analysis and visualization tools 
is already extensively adopted as an interface between Information Technology and Space 
(and suborbital) Science. This is manifested by applications of HEALPix by the following 
projects: CMB experiments Boomerang (deBernardis et al, 2000, Ruhl et al 2003), Archeops 
(Benoit et al, 2003ab), and TopHat, satellite mission WMAP (WMAP2003), the forthcoming 
satellite mission Planck, the Sloan Digital Sky Survey, and others. 

Our involvement with development, distribution, and support of HEALPix since 1997 
would have been impossible without support of a number of institutions and individuals. 
We are indebted to Theoretical Astrophysics Center in Copenhagen, Igor Novikov, and Per 
Rex Christensen; to European Southern Observatory, Peter Quinn, and Kevin Maguire; to 
MPA, Garching, and Simon D.M. White; to Caltech's Observational Cosmology Group, and 
Andrew Lange; to Caltech/IPAC, George Helou, and Ken Ganga; to JPL/Caltech, Center 
for Long Wavelength Astrophysics, and Charles R. Lawrence. We are also grateful for all 
the positive feedback received from numerous HEALPix users worldwide. 
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Table 1: Table of the number of pixels and pixel size accessible to HEALPix . The use of 
32-bit signed integers for the pixel indexing currently restrict the resolution accessible to 
A^ide < 8192. The use of 64-bit signed integers will allow to reach iV sidc = 2 29 . 
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Fig. 1. — Quadrilateral tree pixel numbering. The coarsely pixelised coordinate patch on the 
left consists of four pixels. Two bits suffice to label the pixels. To increase the resolution, 
every pixel splits into 4 daughter pixels shown on the right. These daughters inherit the 
pixel index of their parent (boxed) and acquire two new bits to give the new pixel index. 
Several such curvilinearly mapped coordinate patches (12 in the case of HEALPix , and 6 
in the case of the COBE quad-sphere) are joined at the boundaries to cover the sphere. All 
pixels indices carry a prefix (here omitted for clarity) which identifies which base-resolution 
pixel they belong to. 
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Fig. 2. — cylindrical projection of a number of equal area, iso-latitude tessellations of the 
sphere, which can support a hierarchical tree of further subdivisions of each large base- 
resolution pixel. Six out of possible realizations of such tessellation are shown for several 
values of grid parameters N e , and N^. The HEALPix grid, used for development of the 
corresponding software package, is described by N e — 3, and = 4. 




Fig. 3. — Orthographic projection of the same base-resolution spherical tessellations as those 
shown in Fig. 2. 



Fig. 4. — Orthographic view of the HEALPix partition of the sphere. Overplot of equa- 
tor and meridians illustrates the octahedral symmetry of HEALPix . Light-gray shading 
shows one of the eight (four north, and four south) identical polar base-resolution pixels. 
Dark-gray shading shows one of the four identical equatorial base-resolution pixels. Mov- 
ing clockwise from the upper left panel the grid is hierarchically subdivided with the grid 
resolution parameter equal to N side = 1, 2, 4, 8, and the total number of pixels equal to 
N pix = 12 x N^ ide = 12, 48, 192, 768. All pixel centers are located on N ring = 4 x N side - 1 
rings of constant latitude. Within each panel the areas of all pixels are identical. 
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Fig. 5. — The layout of the HEALPix pixels on the sphere, and demonstration of two possible 
pixel indexetions — one running on iso-latitude rings, and the other arranged hierarchically, 
or in a nested tree fashion. The latter capability is essential for the possible database 
applications of HEALPix. 
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Fig. 6. — Comparison of HEALPix with other tessellations including Quadrilateralized 
Spherical Cube (or QuadCube, used by NASA as data structure for COBE mission prod- 
ucts), icosahedral tessellation of the sphere, and Equidistant Coordinate Partition, or the 
'geographic grid.' The shaded areas illustrate the subsets of all pixels on the sky for which 
the associated Legendre functions have to be computed in order to perform the spheri- 
cal harmonic transforms. This plot demonstrates why the iso-latitude ECP and HEALPix 
points-sets support faster computation of spherical harmonic transforms than the QuadCube, 
the icoshedral grid, and any non iso-latitude construction. 
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Fig. 7. — Illustration of the fundamental difference between computing speeds, which can 
be achieved on iso-latitude and non iso-latitude point-sets. In order to be able to perform 
the necessary computational work in support of the multi-million pixel spherical data sets 
one has to resort to iso-latitude structures of point-sets/sky maps, e.g. HEALPix Moreover, 
the future needs are already fairly clear measurements of the CMB polarization will require 
massively multi-element arrays of detectors, and will produce data sets characterized by a 
great multiplicity (of order of a few thousand) of sky maps. Since there are no computation- 
ally faster methods than those already employed in HEALPix, and global synthesis/analysis 
of a multimillion pixel map consumes about 10 3 s of a standard serial machine CPU time, 
the necessary speed-up will have to be achieved via optimized parallelization of the required 
computing. 



